Telling time with an intrinsically noisy clock 
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Intracellular transmission of information via chemical and transcriptional networks is thwarted 
by a physical limitation: the finite copy number of the constituent chemical species introduces 
unavoidable intrinsic noise. Here we provide a method for solving for the complete probabilistic 
description of intrinsically noisy oscillatory driving. We derive and numerically verify a number of 
simple scaling laws. Unlike in the case of measuring a static quantity, response to an oscillatory 
driving can exhibit a resonant frequency which maximizes information transmission. Further, we 
show that the optimal regulatory design is dependent on the biophysical constraints (i.e., the allowed 
copy number and response time). The resulting phase diagram illustrates under what conditions 
threshold regulation outperforms linear regulation. 



It has long been recognized pQ that the ability to mea- 
sure biochemical quantities, e.g., concentrations, is in- 
trinsically thwarted by the small copy numbers present 
at the scale of the cell. This observation has launched 
considerable experimental investigations as to how high- 
fidelity signal transmission can occur within single cells 
[H [3J , along with an associated literature in mathemati- 
cal and computational techniques for modeling such noisy 
information transmission @H§]. From the perspective of 
biological design - either to understand the mechanisms 
which lead to observed biology or to create synthetic sys- 
tems with desirable properties - these works investigate 
how regulatory elements which comprise biological sys- 
tems function in the presence of intrinsic noise [7] . 

In earlier work we showed how the 'spectral method' 
leads to an efficient and accurate numerical technique, 
which permits optimization to reveal the information- 
optimal design of a transcriptional cascade in the pres- 
ence of intrinsic noise in the statistical steady-state [HI [H] • 
We here turn our attention to the simplest model of the 
dynamic case, illustrated in Fig. [TJa), in which a single 
transcription factor (the 'parent') with copy number n is 
driven by an oscillatory creation rate f(t) = g + a cos uit 
and regulates the expression of a second species (the 
'child') with copy number m; the regulation is modeled 
via the child's creation rate q n . This model captures the 
noisy downstream response to oscillation, e.g., the cell 
cycle, without limiting the results to a particular mech- 
anism for generating oscillations (e.g., via cell division 
[TU] , repressive cycles [3J , or activation-repression circuits 
|llj). We show how the optimal design - i.e., the choice 
of linear- vs. -cooperative and up-vs.-down regulation - is 
determined by the physical demands in terms of allowed 
copy number and response time. Further, while our in- 
tuition from understanding how best to measure static 
signals suggests that slower response time is always more 
accurate [T] , we illustrate how oscillatory driving leads to 
an information-optimal driving frequency, and compute 
how this frequency depends on copy number. 

In the spectral method [H| we exploit the linearity of 
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FIG. 1: (a) A transcription factor (the 'parent') with copy 
number n is driven by an oscillatory creation rate f(t) = 
g + a cos tot and regulates via q n the expression of a second 
species (the 'child') with copy number m. (b, c) Numerical 
confirmations (data points) of analytic expressions (lines) de- 
rived in the small-information limit: Eqns. [9] (circles) and |13| 
top line (squares) and bottom line (triangles) in (b), and Eqn. 
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for both up- (up-triangles) and down-regulation (down- 
triangles) in (c). Parameters are g = no = 1 and c = 0.1 for 
(b), uj = 1 for (c), and a = p = A = 1 and qo = for both, 
yielding small parameters \ui\ < 0.5 (Eqn. EJ, < 0.05 
(linear; Eqn.[l2} and \m \ < 0.184 (threshold; Eqn.[l2l. 



the master equation p nm = - (£„[/(t)] + pC m [q n ])p nm , 
the equation of motion for the joint probability of ob- 
serving n and m copies of the parent and child, respec- 
tively, and expand its solution in terms of the natu- 
ral eigenfunctions of the birth-death process with con- 
stant creation and decay. For a birth-death process 
expressed in terms of an arbitrary creation rate h on 
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species s the positive semidefinite operator C s acts as 
C s [h]p s = sp s - (s + l)]3 s +i + hp s - hp s -i; time is nor- 
malized via the parent decay rate (in these units p is the 
child decay rate). To study dynamics we also Fourier 
transform in harmonics of the driving frequency w, 



In the limit of weak (a <C 1) or fast (w » 1) driving, an 
approximation for /(</>, n) may be obtained by expanding 
in the small parameter \v\\. We first express Eqn. [3] in 
terms of the Fourier transform p(n\<f>) = ~}2 z Pn e ~ lz ^ : 
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where the parent and child eigenfunctions (or 'spectral 

modes') enjoy £ n [g](n\j) = j(n\j) and £ m [q](m\k) = the j = term, i.e. p z 
k(m\k), respectively. Just as in [HI IS], we introduce a 
gauge q to define the basis; the analytic results below are 
independent of this choice. 

The master equation then becomes an algebraic rela- 
tion among the expansion coefficients p z k : 



(7) 

Then we note that for small \v\\, Eqn. p] is dominated by 

z/i|l z l {n\\z\) /\z\\. Since 
this is itself small for z ^ 0, we expand the log in Eqn. [7] 
as log(l + x) — x — x 2 /2 + . . . for small x. The first two 
terms in the log expansion (Sec. [C]) contain the leading- 
order behavior in p z n (proportional to PnP^ 1 = \Pn\ 2 )', 
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where A 3 y = J2n(j\ n )(l ~ 1n)(n\j') . Algorithmically 
we (i) initialize with Gq Q — S z q (set by normalization), 
(ii) exploit the subdiagonality in k, and (iii) for each 
k, exploit the subdiagonality in j; no matrices need be 
inverted. 

Efficient computation of p nm {t) allows optimization 
of the mutual information /(</>, n) between the input 
variable — the phase (j) = u)t of the driving oscillation — 
and the output variable — the copy number of either the 
parent or the child: 
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where p{n\<j>) = p n (t), p{4>) = 1/27T, and p° = 
Ju" d(f) p(n\<j))p((j)) is the time-averaged distribution [T2] , 
Eqn. [3] is integrated numerically during optimization. 

In parallel with the numerical efficiency afforded by 
the spectral method, considerable progress can be made 
analytically. The dynamics of the parent, for exam- 
ple, can be found exactly: the equation for p n (t) (ob- 
tained by summing the master equation over m) is easily 
solved using either the method of characteristics (Sec. 
Al) or spectral decomposition (Sec. A 2). The solu- 
tion is a Poisson distribution with time-dependent mean 
v{t) = Vq + 2\ui \ cos(u;t — 7), where 
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and 7 = tan u). Since the full dynamics are 
known, the Fourier transform coefficients p^ — 
J Q v dip e tz ^ p{n\cf)) / {2-k) are computed by expanding the 
exponential in p{n\(f) and identifying the modes (Sec. 
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where the second to last step uses (n\0) = e 9 g n /n\ and 
(n|l) = (n|0)(n — <?)/<7 fS] to evaluate the sum. Eqn. H 
shows that mutual information asymptotes to the square 
of the amplitude of the oscillation over the mean. Eqn. [9] 
scales like u/° at low frequency and uj~ 2 at high frequency, 
demonstrating that the parent acts as a low-pass filter of 
information; it is numerically verified in Fig. |TJb) . 

Although the child distribution p m (t) is not analyti- 
cally accessible in general, its mean p(t) is exactly cal- 
culable: summing the master equation over both in- 
dices against m and Fourier transforming yields p(t) = 
J2 Z Mz e_iz "S where 
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In the limit of weak regulation, then, (i.e. when q n is 
near constant) we may approximate p m (t) as a Poisson 
distribution with oscillatory mean parameterized by the 
first and second Fourier mode of the exact mean, i.e. 
p(t) ps p,o ± 2\pi \ cos(uit — 6) for up- (down-) regulation, 
where 6 = phasc^^) = tan -1 uj/p + tan -1 u. Under this 
approximation, as in Eqn. [HJ the information between the 
phase of the driving oscillation and the copy number of 
the child is the oscillation amplitude squared over the 
mean, i.e. I((f),m) — \pi\ 2 /po for small \n\\- 

To compare the transmission properties of both non- 
and highly-cooperative regulation, we study both the 
linear function q n = q a + cn and the threshold func- 
tion q n = qg + Ax{n € f&±), respectively, where x is 
a characteristic function equal to 1 when n is in the set 
f2 + = {n > no} (up-regulation) or = {n < uq} 
(down-regulation), and otherwise. In these cases, the 
mean of the child distribution oscillates about the point 
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Here the linear result exploits the fact that the mean 
of the time-averaged parent distribution jp n is g (which 
can be seen from the relationship between distribu- 
tion moments and spectral modes, Sec. ID 1|). In 



the threshold result we define p± = J2nen ± Pn — 

tt±, where ir± = 



^ ngf2 (n|0); the second to last step exploits the result 
EnenJ n \j) = Mno\j - 1) for j > (Sec.p2j and the 
last step takes j = in the small \v\\ limit, 
tude of the oscillation of the child mean is 
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where once more the linear result (top) uses the rela- 



tionship between moments and modes (Sec. D 1 1 and in 



the threshold result (bottom) we define and approximate 
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fi|(n o |0) (Sec. D2). Eqns. 11 and 12 yield the following 



approximations for linear (top) and threshold (bottom) 
regulation: 
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child I(<p, m) is a sharper low-pass filter than the par- 
ent I((j),n), falling off like cj at high frequency in- 
stead of w~ 2 ; it is verified numerically in Fig. [ljb). We 
note that since t, n, and m are not Markov related, i.e. 



p(m\t) ^ ^2 n p(jn\n)p(n\t), Eqn. 13 is not bound by the 
data- processing inequality [13] , and it is possible to have 
I(4>, m) > I((j>, n) (e.g. for linear regulation with uj — > 0, 
<7o = 0, and c > 1), which we have confirmed numerically 
(Sec.|E}. 

Eqn. [13] also offers analytic intuition about the opti- 
mal placement of the parent distribution with respect 
to a threshold regulation function. The derivative of 
Eqn. 13 (bottom) with respect to g vanishes at g* , the 



information-optimal mean of the parent distribution, 
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(recall that dependence on g is contained within (n o |0) 
and 7r±; Eqn. 14 is transcendental and solved iteratively). 
As verified in Fig. |ljc) , Eqn. 14 shows that the par- 
ent distribution is shifted below the threshold for up- 
regulation and above the threshold for down-regulation. 
These shifts account for the ability of up-regulation to 
outperform down-regulation when copy number is highly 
constrained (see Fig. [3] and discussion below) , an effect 
we observed previously [5] when numerically optimiz- 
ing steady-state information between the first and last 
species in a regulatory cascade. 

The above analytic approximations offer guidance dur- 
ing a full numerical optimization of 7(0, m) via the spec- 
tral method. As suggested by Eqn. |13| numerical opti- 
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FIG. 2: (a) At high copy number (N = 20, top curves) the 
optimal information 7* (0, m) exhibits a resonant driving fre- 
quency (point c) for up- (dashed) and down-threshold (dot- 
dashed) regulation, but not for linear (solid) regulation; at low 
copy number (N = 2, bottom curves), there is no resonant fre- 
quency, and slowest (cj — > 0) is best. Panels (b-d) correspond 
to marked points in (a) and show the optimal child distribu- 
tion for down-threshold regulation at phases Q = ut — 8 — 0, 
7r/2, 7r, and 37r/2 [legend in (d) applies to (b-d)]: (b) slow 
driving produces switch-like behavior, with long-lived high- 
(O = 0) and low-copy number (O = n) states and brief inter- 
mediates (0 = 7r/2,37r/2) in between; (c) moderate driving 
produces switch-like behavior with distinguishable intermedi- 
ates, transmitting the most information; and (d) fast driving 
time-averages the parent, and thus the child, distribution. 



mization confirms that m) increases when (i) the am- 
plitude of the driving oscillation is maximal (a — g) and 
(ii) the dynamic range is maximal (<jo = and c — > oo or 
A — > oo). The slope c or discontinuity A, however, is con- 
strained by the average copy number of the child fiQ (Eqn. 
TT| ) . Therefore for a fixed driving frequency and fixed to- 
tal average copy number N = (n) + (m) = g+fio, we opti- 
mize over the single parameter g by setting a — g, q = 0, 
and c = no/ 9 = {N - g)/g or A = — {N — g)/p° ± ; 

additionally we set p = 1 for equal decay rates (as is typ- 
ical when decay rates are dominated by cell division [H] ) . 
For threshold regulation, an optimization over g is done 
at each of a set of values of the (discrete) parameter no, 
and the global optimum is selected. 

At low copy number, optimal information I*((j),m) 
behaves as one might expect from the small-oscillation 
(Eqn. 13 1: it decreases monotonically with fre- 
Ja), bottom curves]. At high copy num- 
ber, I* ((f), m) decreases monotonically with frequency for 
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linear regulation but for threshold regulation exhibits a 
maximum at a resonant frequency [Fig.[2]ja), top curves]. 
Careful examination of the child distribution at different 
phases [Fig. [2](b-d)] (or simply its mean, Sec. [F]) reveals 
the origin of this maximum as follows. As the parent 
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FIG. 3: Phase diagram showing best optimal information 
I*(4>, m) among linear, up-threshold, and down-threshold reg- 
ulation in the space of driving period and copy number. 
Phases are separated by solid lines and marked by sample 
regulation functions (insets). Also shown is (2ir over) res- 
onant frequency u>* [see point c in Fig. [2ja)] as a function 
of copy number for both up- (dashed) and down-regulation 
(dot-dashed). 



N due to its tendency, as discussed above (Eqn. 14 ) and 
in [8] , to require fewer proteins to match the transmission 
across a similarly parameterized down-threshold. 

The low-pass behavior revealed in Eqns. [9] and 
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is 

consistent with our intuition from measuring static quan- 
tities in the presence of intrinsic noise pQ : the longer we 
wait, the more accurate is our estimate. However, in 
the presence of oscillatory driving, we find that threshold 
regulation can lead to an information-optimal frequency, 
and waiting longer is not necessarily the optimal strat- 
egy. Further, we have shown that, at a fixed allowed 
copy number and allowed integration time, one may find 
that a different regulation strategy (linear, threshold up- 
regulation, or threshold down-regulation) is optimal for 
responding to oscillatory driving. Absent from this anal- 
ysis are intriguing questions such as whether the diversity 
of other network topologies observed in nature — includ- 
ing cascades and feedback circuits — are consistent with 
these observations. We anticipate the spectral method 
will continue to be useful in addressing these challenges. 



Appendix A: The parent distribution is a Poisson 
with oscillating mean 



oscillates about the threshold, the child distribution is 
switch-like, with two long-lived switch states centered 
at the threshold's low and high rates, and brief inter- 
mediate states in between. At high copy number, the 
threshold rates are far apart (one is zero and the other 
is large), making the switch states well separated and 
transmitting (slightly more than, due to the intermedi- 
ate states) the strict switch limit [5] of I*(cj>,m) ~ 1 bit. 
For slow oscillations the intermediate states are symmet- 
ric [Fig. |2jb)], but for faster oscillations there is a lag in 
transitioning from one switch state to the other, making 
the intermediate states distinguishable [Fig. |2jc)], and 
transmitting more information about phase. Thus the 
resonant frequency uj* balances the slowness required to 
avoid time-averaging [Fig. [2jd)] with the speed required 
for distinguishable intermediate states. As seen in Fig. 
[3j the onset of a;* occurs above a critical copy number of 
N* ~ 10. 

Fig. |2ja) contains examples in which the most optimal 
regulation function is up-threshold (at low N), down- 
threshold (at high N and high u>) or linear (at high N 
and low lo). The phase diagram (Fig. [3]) shows the re- 
sults of this competition across a range of copy numbers 
N and periods T = 2tt/uj. Linear regulation is best when 
both N and T are large, ultimately surpassing threshold 
regulation's limit of ~1 bit. Down-threshold regulation 
is best at values of T near 2tt juj* because its intermediate 
states are more distinguishable (i.e. have a larger Jensen- 
Shannon divergence) than those of a similarly parameter- 
ized up-threshold. Up-threshold regulation is best at low 



In this study we model two transcription factors - 
the 'parent' with copy number n and the 'child' with 
copy number m — each undergoing a birth-death process, 
with the parent's birth rate an oscillatory function of 
time f(t) = g + a cos cut and the child's birth rate an 
arbitrary function q n of the parent copy number. The 
master equation, 

Pnm — f{t)p n -l,m - f(t)p nrn 
+ (n + l)p n+ i m 

+(m + l)p n,m+l 77lp rxm 

], (Ai) 

describes the time evolution of the joint probability of 
observing n and m copies of the parent and child, re- 
spectively. Here time is normalized via the parent decay 
rate (in these units p is the child decay rate). 

The equation for the parent distribution p n (t) is ob- 
tained by summing the master equation over m: 

Pn = f(t)p n -l - f{t)Pn 

+(n + l)p n +i - np n . (A2) 

Since the parent is not regulated, Eqn. |A2| simply de- 
scribes a one-dimensional birth-death process with time- 
dependent birth rate f(t). The solution can be found, 
regardless of the form of f(i), using either (i) the method 
of characteristics, or (ii) the spectral method; for com- 
pleteness we present both. 



5 



Method of characteristics 



We begin the solution of Eqn. |A2| by denning the gen- 
crating function G(x,t) — ^2 n p n {t)x n P3! over complex 



variable x (writing x 



makes clear that the gener- 



ating function is the Fourier transform in copy number). 
The utility of the generating function is that by summing 
Eqn. |A2[ which describes an infinite set of ordinary dif- 
ferential equations in p n , over n against x n , it becomes a 
single partial differential equation in G, 



G = -(x-l)[d x -f(t)]G. 



(A3) 



The distribution is recovered by inverse transform, 
Pn(t)=d2[G(x,t)} Jnl 

We solve Eqn. |A3| by the method of characteristics, in 
which we demand that x and t lie along a characteristic 
line (x(s),t(s)) parameterized by s, and we seek G(s). 
Expressing Eqn. |A3| as 



(x-l)f(t)G- 



(x-l), 



dG dG 
dt dx 

it is clear that its consistency with the chain rule 

dG _ dG dt dG dx 
ds dt ds dx ds 
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dt 
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(A6) 
(A7) 
(A8) 



dG 

ds v ' dt 

where y = x — 1, the last step of Eqn. |A7| u ses Eqn. |A6| 
and the last step of Eqn. A8 uses Eqns. A6 and A7 We 
integrate Eqn. 



G(y, t) = G exp 



2/o e 
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and recognize that the initial condition Go is an arbitrary 
function of y$, we may therefore expand as Go = J2j c jVo 
for some Cj. Inserting the characteristic equation j/q 
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(Eqn. LA7j, Eqn. A9 becomes 



G(y,t) 
where 
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and the last step in Eqn. |A10| takes to — > — oo to give 
the post-transient behavior, upon which only the j = 
mode survives. Inverse transforming we find 
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(where cq = 1 by normalization), a Poisson distribution 
with time-dependent mean v(t). 

For oscillatory driving f(t) = g + a cos cut, the mean 
evaluates to 



u(t) = g + 



(cos cut + u) sin Lot) (A13) 

1 + co z 



= v a + 2\isi\ cos(uit - 7), 



(A14) 



where = g, 2\vi\ = aj\J\ + oj 2 , and 7 = tan -1 u. Eqn. 
|A14| shows that the parent oscillates about the same point 
and with the same frequency as the driving birth rate, 
but that the oscillation is damped and phase-shifted at 
high frequency. 



2. Spectral method 

Eqn. |A2| can also be solved using the spectral method 
[H [9]. Again we employ the generating function, this 
time expanding in a state space \n) indexed by copy 
number n: \G{t)) = ^2 n Pn(t)\n) . (Projecting onto 
the position space (x\ recovers the previous form with 
(x\G{t)) = G(x,t) and (x\n) = x n .) Summing Eqn. |A2 
against \n) gives 



|G) = -(a+ - l)[fi" - f(t)]\G) = -b+b-(t)\G), (A15) 

where the operators d + and a - raise and lower copy 
number, respectively, i.e. d + \n) = |n+l) and d~\n) — 
n\n — 1) [IMT5] . and we define b + = a + — 1 and b {t) 
11 — f(t). (Note d + f-> x and d~ f> 9^ as is clear from 
Eqns. A15| and|A3"l) 

The spectral method exploits the linearity of the mas- 
ter equation by expanding |G) in the eigenfunctions \j) 
of a birth-death process [8j [9] . Since here the birth rate 
f(t) is time-dependent, we expand as 



(A16) 



where the time-dependent functions (x\j(t)) = (x 



1) 



j i/(t)(x-l) 



are parameterized by the (as yet un- 



A15 



known) function v(t). 
1), the left-hand side of Eqn. 
J2j [Cj \j) + "Cj \j + !>]• Defining !r(t) 
b + f(t) — v{t) such that the \j{t)) are the eigenstates 
of i.e. 



Noting t hat d t {x\j) = v(x\j + 
becomes |G) = 
a - — v(t) = 



b + b-\m) = m)) 



(A17) 



do 



(b + and b raise and lower \j) as a + and a 
respectively), the right-hand side of Eqn. A15 becomes 
Ej Cj + lf(t) - v(t)] \j + 1)}. Therefore project- 



ing onto Eqn. A15 gives the following equation for 
the expansion coefficients Cf. 



-jcj + {f(t) - u(t)] Cj 



(A18) 
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The dynamics are trivial if v = —v + f(i), an equation 
whose solution is Eqn. |All| In this case Eqn. A18 



solved by cj — e - ^* - * ), which becomes Sjq as t — > — oo 
(for post-transient behavior). The probability distribu- 
tion is obtained by inverse transform, 

Pn (t) = (n\G) = £><n|j) = HO) = 

(A19) 

where the last step uses the fact the the zero mode (n|0) 
is a Poisson distribution at the eigenfunction parameter 
(or gauge) v(t) [HI [§]• Eqn. A19 reproduces the result 
from the method of characteristics, Eqn. A12 



Appendix B: Fourier transform of the parent 
distribution 



The Fourier coefficients of a Poisson distribution with 
oscillating mean are here found analytically in terms of 
spectral modes. We begin by representing the distribu- 
tion as 



p(n\4>) = e- u W 



n\ n 



,v{<j>){x-i) 



(Bl) 



where <j> = ojt, and 
v(<f>) = is +2\vi\cos(<t>- 7 ) = ^o+l^l^e^-T). (B2) 



Reordering terms, 



E 



3=0 

x h d * [ eMx ~ 1 H x - i ) j 



X 



E 



=0 
2tt 



£i(j-ey. 

,i(j-2t+z) 



2ir 



(B6) 



The second line in Eqn. |B6| is the derivative represen- 
tation of the spectral mode (n\j) [HI [H|, wit h pa rame- 
ter (or gauge) vq. The fourth line in Eqn. B6 evalu- 



ates to So j- 2i+zi which collapses the I sum nonvanish- 
ingly provided z + j is an even number from to 2j, 
or z G {— j, — j + 2, . . . ,j — 2, j}. This criterion can be 
equivalently be expressed as a condition on j in terms of 
z: j e {\z\, \z\ +2, \z\ +4, . ..,00} = Q, making Eqn.|B£ 

"i\ j (n\j) 



Pn = ^ 



E \(A 



[{j + z)/2]\[{j-z)/2]\ 



(B7) 



Defining j' = (j — \z\)/2 allows the sum to run from to 
00 and yields the simplification [(j + z)/2] \\(j — z)/2]\ = 



j'Kf + \ Z W- f° r au integer z, making Eqn. B7 

|2/+|z| 



Pn 



o lZ ~< 



\- \V\\ 

2-^ oMfV 



j'=0 



(n\2j' + \z\), 



(B8) 



as m Eqn. m 



The Fourier transform will have support only at harmon- 
ics z of the driving frequency, i.e. 



Appendix C: Expansion in the small- information 
limit 



Pn = 



2tt 



Z7T 



2n n! 



(B3) 



= !/ (x-l) 



exp|(a;-l)|z/i|^e ±4W - 7) | . (B4) 



Expanding the exponential and then invoking the bino- 
mial expansion on the ± sum, 



Pn = 



2tv 



dfi iz4> \ an 
2tt n\ x 



o v {x-i) 



E 

3=0 



{x-iy\m\3 



3 .1 

x _ t _. e i(0-7)a-<) <> -i(^-7) 



£=0 



«C7-i)! 



•(B5) 



Here we explicitly expand the log in Eqn. [7] and show 
that the first two terms contribute to the leading-order 
behavior va p n . Eqn. [7] reads 



-iz(j> 



lo s 1 1 + E ^ e " 42 ' 



2tt 



(Cl) 

Since, for small \v\\, Eqn. B8 is dominated by the j' = 
term, i.e. 



Pn 



\A 



\A 



(n\\z\) 



(C2) 



which is small for z ^ 0, we expand the log: 

[■2n 



-iz'cot 




(C3) 
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or, explicitly writing out the I sums, 



Appendix D: Useful properties of the spectral modes 



(C4) 



Reordering terms as 



1 = E 



(-1) 



e+i 



T — 

L^i Cr.0 



(p°nY 



1=1 

Zl^O z 2 #0 z 

x / ^(p i ( z+Zl+Z2+ ... +ze ) < t > 

Jo 2tt 



(C5) 



and employing Jg 2 ^ d</> e %s< ^ — 27T(5 s o for integer s allows 
us to collapse the z sum, 



1 = E 



E- 



£=1 
X 



E P« E Pn ' ' ' E Pn 

(zi+z 2 H hZ£) 



(C6) 



Now, since p n oc for small |^i| (Eqn. |C2]), I(<j>,ri) 

will be dominated by the leading-order term in \z\. Writ- 
ing out the first few terms in Eqn. |C6| explicitly, 



1 = E^oE^ 



2 Z-^ I'll "; 2 /_^Pn /_^Pnh>n 

n Wn) x^O y^O 

+ 3E7^EKEKE^ ( " + ^ 
+ ..., (C7) 

we see that the leading-order term is proportional to 
PnPn 1 = |Pn| 2 an d nas contributions from both the first 
(z = ±1) and the second ((x,y) = (±1, =pl)) term in the 
log expansion. To leading order, then, 



1 + 1-^(1 + 1) 



V- bnl 2 V- \pI\ 2 ,p„. 

2^—o- = 2^—o-> ( L8 ) 

r, Pn „ Pn 



as in Eqn. [8] 



1. Relating spectral modes to distribution 
moments 

The eigenfunctions of a birth-death process, or 'spec- 
tral modes,' define a complete basis in which any distri- 
bution p n can be expanded, i.e. 



P 



» = E^'Hj) 



(Dl) 



(where the expansion coefficients are computed by inverse 
transform: Cj — J2 n Pn(j\ n ))- The goal of this section is 
to relate the spectral modes to the distribution moments 

The £th moment is (n e ) = ^ ■ CjT)j, where 



^£^< n |i>=E'< 



dx (x - \ye a{x -^ 



2ni 



-n+l 



(D2) 



The last step uses the integral representation of (n\j) over 
complex variable x [9]. Isolating the n-dependence, we 
may write 

V^i^ix-iye^-Se, (D3) 
J J 2m x 

where for w = 

5, = V f»V* = V «;» = {wd w f — !— , (D4) 

z — / z — / 1 — w 

n n 

and the last step sums the geometric series. Computing 
the first few derivatives in Eqn. |D4| reveals the pattern 



in terms of the Eulerian numbers 

7/. 

+ 1 

v=0 

making Eqn. |D3| 



dx 



2tt« (a; - l^-i+ie-sO" 1 ) ' 



(D7) 



The integral is recognized as a representation of the con- 
jugate mode (j'\n') with f = £ — j > 0, n' = u > 0, 
and parameter (or gauge) — g [§]• Modes and conjugate 
modes can be evaluated either recursively using selection 
rules or explicitly using Cauchy's theorem Conjugate 
modes (j'\n!) with gauge g' are j'th order polynomials in 
n'; the first few are 

(j' = 0\n') = 1, (D8) 
(j' = lK) = n'-g', (D9) 

0' = 2K) = in' 2 -r 9 ' + i)n' + ^. (D10) 
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Thus the first few m are (Eqn. D7 1 

$>li) = 7f° = *Jo, 

n 



(Dll) 
(D12) 



= 5 («? + l)^o + (2.g + l)^i + 26 j2 , (D13) 

and in general the £th moment is calculated (n £ ) — 
£j r j'ir 



Eqns. Dll and D12 are useful in describing the mean 
properties of a child species regulated by the linear func- 
tion q n = qo + cn. As described in Eqn. [TTJ the mean of 
the child distribution oscillates about the point 



Mo = £ q nP Q n = £(go + cn) (n ->./ : • ( D I 4 1 



2.; 



where the last step uses the Fourier transform of the par- 
ent distribution, Eqn. |B8[ Evaluating the n sum using 
Eqns. [DTT1 and |Di~2| gives 



Mo 



E 



l 2 J 



(JO 2 



= qo + eg, 



[qo$2j,o + c{gS 2 jfi + <%i)] (D15) 

(D16) 



(since only the first and second Kronecker delta have sup- 
port for integer j), as in Eqn. 11 (top). 

Similarly, as described in Eqn.|12| the amplitude of the 
oscillation of the child mean is 



iMil 



X^2{q +cn)Y,- 



|2j + l 



(n\2j + l). (D17) 



Again employing Eqns. Dll and |D12} 



1 



Mil 



|2j + l 



v/l + ^/p) 2 + 
x + c(3<5 2: ,+i,o + <%+i,i)] (D18) 

(since only the last Kronecker delta has support for inte- 
ger j), as in Eqn. 12 (top). 



2. Sums of differences 

The zero mode (i.e. the steady-state solution) of the 
birth-death process with birth rate (or gauge) g is the 
Poisson distribution, 



(n\j=0) =e-«- 



(D20) 



(recall that time is normalized by the decay rate). Each 
higher mode is related to the previous mode by discrete 
derivative [5]: 



(n\j + 1) = -d n (n\j) = (n - l\j) - (n\j). 



(D21) 



This property is especially useful when performing a fi- 
nite sum, since only the boundary terms survive: 



j>|i + l) = <a-l|i)-<&|j>. 



(D22) 



We make use of Eqn. |D22| in arriving at Eqns. [TT] 
(bottom) and 12 (bottom), describing the mean prop- 



erties of a child species regulated by the threshold func- 
tion q n = qo + A%(n G f2±). Here x is a character- 
istic function equal to 1 when n is in the set J7 + = 
{n > n } = {n + l,...,oo} (for up-regulation) or 
f2_ = {n < no} = {0, ...,no} (for down-regulation), 
and otherwise. In Eqn. |TT| (bottom), 



E QnPn 



■2.1 



(j!) 2 



'7o 



l 2 J 



j KJ '> n£(l± 



(D23) 
(n|2j),(D24) 
(D25) 



where the go term collapses as in Eqn. |D16| The 
j = term reduces explicitly to Air±, where n± = 
SnGf2+ For an y f > we may employ Eqn 



D22 



E ("IJ") = (D26) 

n— no + 1 

5>|j') = -(no|i'-l), (D27) 

n=0 

(since the boundary terms at n = — 1 and 00 vanish), 



making Eqn. D25 



Mo = go + A 



3>0 



(j!) 2 ' 



-±±EfeW(«o|2. ? -l> 



as in Eqn. [TT] (bottom). In Eqn. 12 (bottom), 



Mil 



1 

TiTWpF 



|2j + l 



x ?STIji ( " |!i + 1 » 



(D28) 



(D29) 



(D30) 



A 



E|S 



V / 1+W?Y J!(J + I)! n 
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parent information: 
l((|),n) = 0.368 bits 
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child information: 
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O 8 = 3n/2 



10 15 
child copy number, m 



•a § @ 8 1 a a a a a 
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FIG. 4: The child can tell better time than the parent. Parent 
(a) and child (b) distributions are plotted at phases V = cut—"/ 
and O = Ldt—9 equal to 0, it/2, n, and 3ir/2, where 7 and 6 are 
the phase shifts of the parent and child means, respectively. 
Regulation is linear and steep (slope c = 3); other parameters 
are g = a = p = 1 and uj = go = 0. 



where the qo term vanishes as in Eqn. |D19| Again using 
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Eqns. D26 and D27 



2 4 6 

driving phase, <J> = cot 



l"i 



2j+l 



i!(i + i)! 



K|2j>, 



(D32) 



as in Eqn. 12 (bottom), where for the down-threshold 
the negative sign is absorbed into the definition /i(t) « 
^0 ± 2|^ii| cos(ujt — 9). 



FIG. 5: Information-optimal child distribution mean /i(i) = 
~}2 m Pm(t) (solid) as a function of phase c/> = u)t for threshold 
down-regulation and fixed copy number N = 20 and driving 
frequency w = (a) 0.10, (b) 0.21, (c) 0.46, (d) 0.91, (e) 1.91, 
and (f) 3.98. Panel (c) corresponds to resonant frequency 
ui* at which optimal information I*(cj>, m) is greatest. For 
comparison, dashed lines show square waves centered at no 
with amplitude 2|^i| and phase shift 6 = phase(/ii). 



3. Differentiation with respect to gauge 



rTil which uses 



The derivative of spectral mode (n\j) with respect to 
its gauge g is 



d g (n\j)=(n\j + l), 



(D33) 



a property that can be seen most readily from the integral 
representation of (n\j) (Eqn. D2 1. This property is useful 



in expediting derivative calculations, for example Eqn. 



d g (n o \0) = (no|l> = (no|0) 



no 



(D34) 



d g ^± = d g J2 H°)= E Hl) = ±(n |0),(D35) 



where the last step in Eqn. |D34| uses Eqn. D21[ and the 
last step in Eqn. |D35|uscs Eqns. ID26I and ID27 
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Appendix E: The child can tell better time than the 
parent 



Appendix F: Threshold regulation can exhibit a 
resonant frequency 



Although the master equation (Eqn. Al I is Marko- 



vian in time (i.e. the probability of making a transition 
at time t is independent of previous transitions), it is 
not explicitly Markovian in the variables t, n, and m: 
p(m\t) p(m\n)p(n\t) . As such, information trans- 

mission is not bound by the data-processing inequality 
1 1 3 j . and it is possible for the child to transmit more 
information than the parent about the driving phase, 
I((f>,m) > 7(0, n). This possibility is explicitly appar- 
ent in the small-oscillation limit with linear regulation, 
Eqn. [13| (top), 



(q + cg)[l + (u]/ P y 



-J{(j>,ri), 



(El) 



for example if uj — > 0, <?o = 0, and c > 1. However, be- 
cause these are the very parameter settings that strain 
the approximations under which Eqn. El is derived (i.e. 



weak or fast oscillation and near-constant regulation), it 
is useful to also demonstrate numerically via the spectral 
method a case in which I(<f),m) > I(cf>,n). Fig. [1] shows 
clearly that if the regulation is sufficiently steep, the os- 
cillation is sufficiently amplified and the child tells better 
time than the parent does, i.e. I(4>,m) > /(</>, n). 



When the regulation function is a threshold and the 
total copy number N is sufficiently high, the optimal in- 
formation /*(</>, m) exhibits a maximum at a resonant 
frequency cu* . In Fig. [2] it is shown that u>* is the fre- 
quency at which the driving oscillation is slow enough for 
the output to be switch-like (and avoid time-averaging), 
but fast enough for the brief states in between the switch 
states to be distinguishable from each other. 



Fig. [5] here plots the mean fi(t) of the child distribu- 
tion against phase for a range of driving frequencies. At 
low frequency [Fig. [5ja)], the output is switch-like, and 
fi(t) is well approximated by a square wave. At the res- 
onant frequency [Fig.^c)], the output is still switch-like 
but no longer symmetric in time, the asymmetry arising 
from a lag in transitioning from one switch state to the 
other (the lag allows the transitions to be distinguished 
from each other, maximizing transmission of information 
about phase). At high frequency, [Fig. [5^f)] the driving 
is faster than the parent decay rate, and both parent and 
child distributions are time-averaged. 
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